Stability of two-dimensional spatial solitons in nonlocal nonlinear media 
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We discuss existence and stability of two-dimensional solitons in media with spatially nonlocal 
nonlinear response. We show that such systems, which include thermal nonlinearity and dipolar Bose 
Einstein condensates, may support a variety of stationary localized structures - including rotating 
spatial solitons. We also demonstrate that the stability of these structures critically depends on the 
spatial profile of the nonlocal response function. 
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I. INTRODUCTION 

A spatial optical soliton is a beam which propagates in 
nonlinear medium without changing its structure. For- 
mation of solitons is a result of a balance between diffrac- 
tion caused by the finite size of the wave and the non- 
linear changes of the refractive index of the medium in- 
duced by the wave itself When this balance can be 
maintained dynamically the soliton exists as a robust ob- 
ject withstanding even strong perturbation. It appears 
that solitons, or solitary waves, are ubiquitous in nature 
and have been identified in many other systems such as 
plasma, matter waves or classical field theory. 

Solitons have been typically considered in the context 
of the so called local nonlinear media. In such media the 
refractive index change induced by an optical beam in a 
particular point depends solely on the beam intensity in 
this very point. However, in many optical systems the 
nonlinear response of the medium is actually a spatially 
nonlocal function of the wave intensity. This means that 
index change in particular point depends on the light 
intensity in a certain neighborhood of this point. This 
occurs, for instance, when the nonlinearity is associated 
with some sort of transport processes such as heat con- 
duction in media with thermal response ^0,0) diffusion 
of charge carriers J] or atoms or molecules in atomic 
vapors lii • It is also the case in systems exhibiting 
a long-range interaction of constituent molecules or par- 
ticles such as in nematic liquid crystals hOl LlL ^|| 
or dipolar Bose Einstein condensates 0, llSl Ha llTlllSj. 
Nonlocality is thus a feature of a large number of non- 
linear systems leading to novel phenomena of a generic 
nature. For instance, it may promote modulational in- 
stability in self-defocusing media (I&. i2(i[2ll l23|. as well 
as suppress wave collapse of multidimensional beams in 
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self- focusing media [Sl|23,|2J|. Nonlocal nonlinearity may 
even represent parametric wave mixing, both in spatial 
j25| and spatio-temporal domain 26] where it describes 
formation of the so called X-waves. Furthermore, non- 
locality significantly affects soliton interaction leading to 
formation of bound state of otherwise repelling bright or 
dark solitons "27, "2^ |2^ • It has been also shown that 
nonlocal media may support formation of stable complex 
localized structures. They include multihump [sil Is^ l 
and vortex ring solitons [33, |3J| . 

The robustness of nonlocal solitons has been attributed 
to the fact that by its nature, nonlocality acts as some 
sort of spatial averaging of the nonlinear response of the 
medium. Therefore perturbations, which normally would 
grow quickly in a local medium, are being averaged out 
by nonlocality ensuring a greater stability domain of soli- 
tons 33, 35] . However, it turns out that such an intuitive 
physical picture of nonlocality-mediated soliton stabiliza- 
tion which has been earlier confirmed in studies of mod- 
ulational instability and beam collapse, has only limited 
validity. Recent works by Yakimenko [s^], Xu [3^ and 
Mihalache |33] demonstrated that even a high degree of 
nonlocality may not guarantee existence of stable high 
order soliton structures. 

In this work we investigate the propagation of finite 
beams in nonlocal nonlinear media. We demonstrate that 
while nonlocality does stabilize solitons, the stability do- 
main is strongly affected by the actual form of the spatial 
nonlocality. 

This paper is organized as follows: In Section ^ we 
introduce the mathematical models under consideration 
and the ansatz for the solutions we are interested in. Sec- 
tion ^O] deals with spatial optical solitons in media with 
a thermal nonlinearity. A remarkably robust new type 
of rotational soliton is presented. In Section Hvl we treat 
the more complicated case of two-dimensional solitons 
in dipolar Bose Einstein condensates (BEC). Due to the 
mixture of local and nonlocal type of nonlinearity a vari- 
ety of solutions can be stabilized. Finally, in Sect ion Ivl we 
discuss the observed dynamics of the soliton structures. 
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II. MODEL EQUATIONS 

We consider physical systems governed by the two- 
dimensional nonlocal nonlinear Schrodinger (NLS) equa- 
tion 

where 9 represents the spatially nonlocal nonlinear re- 
sponse of the medium. Its form depends on the details 
of a particular physical system. In the following we will 
consider three nonlocal models. 

The first one is rather unphysical but extremely in- 
structive, the so called Gaussian model of nonlocality. In 
this model 9 is given as 

^=^jj (2) 

In spite of the fact that there is no known physical sys- 
tem which would be described by a Gaussian response, 
this model has served as a phenomenological example of 
a nonlocal medium, enabling, thanks to its form, an an- 
alytical treatment of the ensuing wave dynamics. 

The second model, referred to as thermal nonlinearity, 
describes, for instance, the effect of plasma heating on the 
propagation of electromagnetic waves ^ . In this case 9 
is governed by the following diffusion- type equation 0, H 

/ 52 ^2 \ 

which is valid for the typical spatial diffusion scale large 
compared to the operating wavelength • Solving 

formally Eq. (j^J in the Fourier space yields 

0=^ JJ M\r~r^\)mr',z)fd'r^, (4) 

where ^ is the modified Bessel function of the second 
kind and r = xex + ySy denotes the transverse coordi- 
nates. 

The third system considered here is the model of a 
dipolar Bose Einstein condensate where the nonlocal 
character of the interatomic potential is due to a long- 
range interaction of dipoles [22|- Such condensate has 
been realized recently in experiments with Chromium 
atoms which exhibit strong magnetic dipole moment [l5j | . 
Considering only two transverse dimensions and time 
(which plays a role analogous to that of the propagation 
variable z) one arrives at the following formula [18i | 

^ = "IV'I' + ^ JJ^{\r-r^\)mr',z)fd^r^, (5) 

and 

m = JJ W^|^f^e^(^-+^«^)dfc.dA„ (6) 




Transverse Coordinate if 

FIG. 1: Nonlocal response functions for thermal nonlinearities 
(.So) and dipolar EEC's (JH). For comparison, a Gaussian 
response function is shown, too. 



with k = y fc^ + and erfc being the complementary 

error function 0, |43| • Interestingly, this model contains 
both, local and nonlocal components. The parameter a 
controls the sign as well as the relative strength of the 
local component of nonlinearity. This model is somewhat 
controversial and its validity has been recently questioned 
41]. However this issue is beyond the scope of this paper 
and we consider Eq.® here in its own right, as another 
example of nonlocal nonlinearity. 

Note that the above three models are represented in 
dimensionless form. This means that only the spatial 
extent of the actual solution determines whether one 
operates in a "local" (|'0P varies slowly over the trans- 
verse lengths of the order unity) or "strongly nonlocal" 
(IV'P changes fast) regime. The "conventional" degree of 
nonlocality, the actual width of the response function, is 
scaled out. 

It is known that the thermal model [Eq. ^] permits, 
apart from the ground state, stable single-charge vortex 
solutions only, whereas multi-charge vortices appear to 
be unstable (34] . In contrast, when a Gaussian response 
function is used [Eq.Q], also stable multi-charged vor- 
tices may exist jSSj. Figure ^ illustrates the nonlocal re- 
sponse functions for the two physically relevant systems 
introduced above, and, as a comparison the (unphysi- 
cal) Gaussian response. It is obvious that the thermal 
and condensate response functions exhibit qualitatively 
the same spatial character. Therefore, one can expect a 
similar behavior of solitary solutions for beams in media 
with thermal nonlinearities and matter waves in dipolar 
EEC's, especially with respect to their stability. 

In the following analysis we will be interested in two 
classes of nonlinear localized solutions: classical solitons 
with stationary intensity and phase distribution, and a 
new recently introduced type of rotating solution, the so 
called " azimuthons" • In case of the standard solitons 
we are looking for stationary solutions in the form 

V'(r,(/.,z) = C/(r,0)e^^^ (7) 

where A is the propagation constant. Here we changed 
from Cartesian to cylindrical coordinates for technical 



convenience. In the highly nonlocal limit the soliton pro- 
file |C/p(r, (/)) is expected to be very narrow compared to 
a width of the response function (which is a fixed quan- 
tity in our dimensionless system). Under this condition, 
for the Gaussian response the nonlinear term Eq. (0) can 
be represented in simplified form |24j 



= — Pocxp(-rV2), 



(8) 



where Pq — J \U\'^(Pfis the total power. Now the nonlin- 
ear Schrodinger equation Eq. JQ) becomes linear and lo- 
cal. It describes the evolution of an optical beam trapped 
in an effective waveguide structure ("potential") with the 
profile given by the Gaussian nonlocal response function. 
This highly nonlocal limit was first explored by Snyder 
and Mitchell in the context of the "accessible solitons" 
141^ . In our scaling the strongly nonlocal limit corre- 
sponds to high power Pq —^ oo. Large power Pq means 
now a deep "potential" —9 in Eq. ^ and stronger con- 
finement of modes (solitons) in the vicinity of r = 0. We 
exploit this accessibility character of the solitons in our 
numerical scheme (see Appendix IXj) . 

The second class of solutions we deal with, the az- 
imuthons, are a straightforward generalization of the 
ansatz {Tj). They represent spatially rotating structures 
and hence involve an additional parameter, the rotation 
frequency il 



ip^r, (j),z) — U (r, (j) — r2z)e' 



(9) 



For = 0, azimuthons become ordinary (nonrotating) 
solitons. For example, the most simple family of az- 
imuthons is the one connecting the dipole soliton with 
the single charged vortex soliton (for fixed propagation 
constant A). The single charged vortex consists of two 
dipole-shaped structures in real and imaginary part of U 
with equal amplitude. If these two amplitudes are dif- 
ferent we have a rotating azimuthon, if one of the ampli- 
tudes is zero we have the dipole soliton. In the following 
we will refer to this amplitude ratio in terms of the mod- 
ulation depth 



jmaxReC/ — maxlm[/| 
max \U\ 



(10) 



In case of a Gaussian nonlocal response function, for 
instance, the azimuthons arise naturally again via the 
fact that — Pq exp(— r^/2)/27r for high powers. Since 
Eq. becomes linear, solitons will converge to the lin- 
ear modes of the corresponding "potential" —9. In this 
strongly nonlocal limit azimuthons, as introduced above, 
will converge to two degenerate dipole modes with the 
phase difference of 7r/2 and unequal amplitudes. As 
pointed out in |42l |. the rotation frequency f2 is deter- 
mined by the amplitude ratio of the two degenerate 
modes (modulation depth) and, of course, the propaga- 
tion constant A. 
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FIG. 2: Dynamics of tlie vortex states in tliermal model: (a) 
unstable propagation at Pq = 200, profiles are shown in a xy 
box of about 5x5; (b) unstable propagation at Pq = 500, 
profiles are shown in a xy box of about 3x3; (c) stable at 
Po = 2000, profiles are shown in a xy box of about 1.5 x 1.5. 



III. BEAM PROPAGATION IN MEDIA WITH A 
THERMAL NONLINEARITY 



Let us first have a look at Eqs. ^ and (0J. It is known 
that in the local classical NLS, no stationary state exists. 
Finite beams either diffract if their power is too low or 
experience catastrophic collapse if their power exceeds a 
certain threshold. However, it has been already shown 
that in case of nonlocal NLS the catastrophic collapse is 
arrested and fundamental-type solitons can be stabilized. 
Moreover, it turns out that for sufficiently high degree of 
nonlocality (in our scaling this is equivalent to sufficiently 
high power), even the single charged vortex is reported 
to become stable [S^. In fact, it is so far the only known 
stable stationary solution in this system, apart from the 
fundamental, single peak soliton. Indeed, Fig. [21 shows 
that for power Pq = 2000 the single charged vortex-state 
is stable, while for Pq = 200 and Pq = 500 it decays 
after a certain propagation distance. Of course, strictly 
speaking we cannot claim stability from numerical prop- 
agation experiments, since we can propagate beams over 
finite distances only. However, from our results we can 
at least infer that the potential growth rates of unstable 
internal modes are very small. For instance, the vortex 
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FIG. 3: Dynamics of azimuthons in thermal model of non- 
linearity. Modulation depth n — 0.5: (a) stable propagation 
at Po = 500, rotating with S7 ~ 4.9, profiles are shown in a 
xy box of about 3x3; (b) stable propagation at Po — 2000, 
rotating with Q 11, profiles are shown in a xy box of about 
i.5 X i.5. 
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FIG. 4: Dynamics of localized structures in a BEG nonlin- 
ear model [Eq.0]. Power Po = 1000, a — 0: (a) unstable 
vortex-state; (b) stable azimuthon, rotating with fl ~ 5.1, 
modulation-depth n = 0.5. All profiles are shown in a xy box 
of about 4x4 



state in Fig.|2Ic) was stable over more than 10^ numerical 
2;-steps. The small periodic oscillations visible in some of 
the curves of numerical figures are due to the excitation 
of stable internal modes of the solution by perturbations 
of the initial data. In fact, since we use approximate so- 
lutions computed by the method described in Appendix 
IXI we always deal with perturbed initial data. 

Both vortices with Pq — 200 and Po = 500 are unsta- 
ble, but their decay dynamics are quite different. The 
one with lower power [Fig. [2Ia)] breaks up into two 
ground state solitons, which move away from the cen- 
ter once formed which is a typical instability scenario 
for local NLS. In contrast, for higher powers (or, con- 
versely - stronger nonlocality) the vortex initially breaks- 
up and later fuses to a single ground state solution emit- 
ting remnants of its original structure [seel^^b)]. This 
behavior can be explained by the nonlocality-induced 
broad single waveguide (or attractive potential) which 
prevents fragments from leaving (unlike the break-up in 
local medium). 

The decay behavior of the Po — 500 vortex indicates an 
effective power flux towards the center. Hence, one might 
expect that a rotation of the object could somehow bal- 
ance this contraction. As mentioned in Sec. ^ we are 
not looking at classical solitons only, but also at rotating 
azimuthons. At each power level Pq we expect a fam- 
ily of azimuthons, allowing a continuous transition from 
the single charged vortex (modulation-depth n = 0) to 
the flat-phased dipole-state (modulation-depth n = 1). 
It turns out that the idea that the rotation of the az- 
imuthon might balance the contractive forces destroying 
the vortex is indeed correct [s^ . Fig. Olb) shows that an 
azimuthon can be stable at power levels, where normally 



the vortex decays collapsing into a fundamental soliton 
(Po = 500). This shows that azimuthons are more robust 
than vortices in this system. Hence, the azimuthons may 
be easier to realize than vortices in future attempts of ex- 
perimental observation of higher order nonlinear beams 
in media with thermal nonlinearity. 

In our extensive numerical simulations of the thermal 
model we were not able to observe any other stable, non- 
rotating dipole-states. Even at higher powers (stronger 
nonlocality) the simple dipole-state remained unstable. 
Moreover, we did not observe any decrease of the ob- 
served growth rates with increasing powers (shift of the 
break-up point to larger values of z). Also higher or- 
der states like quadrupole etc. were found to be always 
unstable in this system. 



IV. TWO-DIMENSIONAL DIPOLAR 
BOSE-EINSTEIN CONDENSATES 

Compared to the previous system, the system of Eqs. 
and jSJ is more complicated. A free parameter a is 
involved, which controls both, the sign and the relative 
strength of an additional local contribution to the non- 
linear response. Let us first consider the case a = 0. It 
turns out that in this case the system behaves qualita- 
tively like the thermal model discussed in the preceding 
section. The azimuthons are found to be the most ro- 
bust objects (apart from the stable ground state), and 
we observe the same stabilizing-by-rotation mechanism 
as above. The vortex at Pq — 1000 shrinks to form a 
single ground state solution [Fig. 2Ja)], whereas the ro- 
tating azimuthon at the same power-level is stable due to 
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FIG. 5: Stable vortex-states in the BEG model: (a) stable 
propagation at Po = 5000, a = 0, profiles are shown in a 
xy box of about 2x2; (b) stable propagation at Po = 1000, 
a = —0.05, profiles are shown in a box of about 5x5. 



its rotation [Fig.^b)]. When we increase the power suf- 
ficiently the vortex-state also becomes stable [Fig. |3Ia)], 
as expected. 

The local part of the nonlinearity [Eq.(|SJ|], can be ei- 
ther of focusing or defocusing nature, depending on the 
sign of the coefficient a. If a is positive, stable solu- 
tions become unstable for a large enough. E.g., if we fix 
the power Pq = 1000, the azimuthon with modulation- 
depth n — 0.5 [stable for a — 0, see Fig. 01b)] turns 
out to be unstable already at a = 0.01. On the other 
hand, if a is negative, the otherwise unstable solutions 
can be stabilized. Fig.^Ib) shows the stable vortex-state 
at Po = 1000, a — —0.05, but the corresponding vortex 
state at a = is unstable [see Fig.QJa)]. Moreover, while 
playing with the parameter a we might even stabilize the 
dipole-state, as shown in Fig. (SJb). This is remarkable, 
since the dipole seems to be always unstable if a = 
[Fig. Eta)], even for higher powers, and also in the sys- 
tem considered in Sec. IIIII However, the parameter a 
needs to be tuned carefully. If it is too large in absolute 
value, the stabilization mechanism fails [see Fig.|H|[c)]. 

In a last example we show that even higher order soli- 
tons which are not related to the single-charged vortex 
(like the dipole) can be stabilized. The double charge 
vortex turns out to be stable with a small defocusing 
local contribution {a — —0.03), as revealed by Fig. [7| 
The systematic studies of stability of other higher order 
solitons is beyond the scope of this paper. However, it 
seems that systems featuring a strong focusing nonlocal 
nonlinearity combined with a small local contribution of 
defocusing nature are good candidates for observation of 
stable higher order solitons. 
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FIG. 6: Dipole states of the BEG nonlocal model at Pq = 
1000: (a) unstable propagation, a — 0, profiles are shown in 
a xy box of about 4x4; (b) stable propagation, a = —0.03, 
profiles are shown in a xy box of about 5x5. (c) unstable 
propagation, a = —0.05, profiles are shown in a box of about 
5x5. 
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FIG. 7: Dynamics of double charge vortex states in a BEG 
nonlocal model with Po = 2000: (a) unstable propagation, 
a — 0, profiles are shown in a xy box of about 4x4; (b) 
stable propagation, a = —0.03, profiles are shown in a xy box 
of about 4x4. 
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V. DISCUSSION 

The examples considered above clearly demonstrate 
that the stability of higher order nonlinear modes cru- 
cially depend on the shape of the nonlocal response func- 
tion. Briedis et al. \^?^ suggested a very illustrative ex- 
planation why for a Gaussian response function higher or- 
der solitons can be stable. In the strongly nonlocal limit 
(high powers), the soliton shape U becomes very nar- 
row compared to the Gaussian response function. Hence, 
in this hmit Eq. |(2Jl simphfies to 6* = Pq exp(— r^/2)/27r 
where Pq = J \ U\'^cfif, and Eq. Q becomes linear and 
local. Hence we expect higher order states become stable 
at high enough powers in the Gaussian system. We may 
get at least some idea why in our physical systems [Eqs. 
Q and ISJ] higher order solitons, e.g. multi-charged vor- 
tices, are never stable: If we have a look at Fig. ^ it is 
obvious that the argument of Briedis et al. [33| breaks 
down for the response functions .^o and JH. The singu- 
larity of these functions at r = does not permit the 
simplification of the convolution integral like in the case 
for the Gaussian response function. No matter how nar- 
row the soliton shape U may be, it will always feel the 
singularity of the response functions .^o and 91. 

Even the role of the local nonlinearity (parameter a) 
in the dipolar BEG [Eq. can be understood in this 
context. If a is positive, the focusing local contribution 
makes the singularity of the response-function [Eq. |(5j] 
at r = "worse", and therefore tends to destabilize non- 
linear solutions. On the other hand, for negative a the 
defocusing local contribution may balance the destabiliz- 
ing effect of the nonlocal part of the response-function. 

A similar observation was reported recently by Xu et 
al. for the one-dimensional nonlocal NLS equation. 
In the ID case, a thermal nonlinearity like Eq. Q leads 
to an exponential response function exp(— 2|x|), where 
X is the transversal coordinate. In the strongly nonlo- 
cal regime for this exponential response only fundamen- 
tal, dipole-, triple- and quadrupole-states are reported to 
be stable, whereas for the Gaussian response all higher 
order states turn out to become stable if the nonlocal- 
ity is sufficiently high. Again, we see that due to the 
non-continuous first derivative of exp(— 2|a;|) the Briedis 
et al. argument for stability does not hold. For com- 
pleteness we investigated numerically several (unphysi- 
cal) response functions in ID case. It seems that smooth 
response functions like exp(— x^), sech(a:), 1/[1 + (x)'^] 
etc. permit stable multipole-states of high order in the 
strongly nonlocal limit. In contrast, response functions 
with an apex at a; = like exp(— 2|a;|), 1/(1 -I- 2|x|)^, 
or (n -I- 1)(1 — ■\/[x[)/2 with < 1, ri > 1 allow only 
stable monopole-, dipole-, triple- and quadrupole-states 
in this regime. If we choose a response function with a 
singularity at x = 0, hke {l/^/\x\ — l)/2, at least also the 
quadrupole-state becomes unstable in the highly nonlo- 
cal regime. However, we can state that both in 2D and 
ID, once the Briedis et al. argument breaks down, even 
a strong nonlocality may not stabilize solitons of suffi- 



ciently high order. 

VI. CONCLUSIONS 

In conclusion we discussed the propagation of 
two-dimensional solitons in nonlocal nonlinear media. 
We considered two physical relevant systems, optical 
beams in media with a thermal nonlinearity and two- 
dimensional dipolar Bose-Einstein condensates and com- 
pared them with the pure Gaussian nonlinear response 
which seems to support variety of stable solitons of high 
order. We demonstrated that while nonlocality does 
stabilize solitons, the stability domain is strongly af- 
fected by the actual form of the spatial nonlocality. We 
showed that both physically relevant systems support 
stable propagation of a new type of rotating solitons. 
In fact, our results suggest that these rotating structures 
are the most stable higher order nonlinear states in these 
systems, and might be even easier to realize experimen- 
tally than a single charge vortex or dipole solitons. Their 
appropriate input intensity and phase spatial profiles can 
be easily prepared using, for instance, spatial light mod- 
ulator. We also showed that in a particular case of the 
BEG system, the mixture of local and nonlocal response 
stabilizes some higher order localized states in certain 
parameter regions. 
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APPENDIX A: NUMERICAL ALGORITHM 

In order to find stationary localized structures sup- 
ported by the nonlocal models Eqs. 13 (|2J) one must 
solve the 2-dimensional nonlinear Schrodinger equation 
with corresponding nonlinear nonlocal term. In general, 
it is a quite difficult task to compute stationary solutions 
in 2D, especially when one cannot exploit certain sym- 
metries, like, e.g., cylindrical symmetry as in the case 
of vortices. Here, we present a very fast and easy way 
to compute at least approximate solutions which relies 
on, the so-called, strongly nonlocal limit which we will 
illustrate using the Gaussian model. 

In the strongly nonlocal limit finding the exact sta- 
tionary solution of the propagation equation reduces 
to the eigenvalue problem Av = Xv, where vector v 
represents discretized function U{r) and A is a sym- 
metric band matrix. These facts can be exploited to 
compute selected eigenvalues and eigenvectors iteratively 
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FIG. 8: Approximate soliton solutions of Eqs. Q and ^ 
for Po = 1000: (a) dipole state Uo computed from 9o — 
Po exp(— r'^/2)/27r; (b) dipole state C/g after six iterations 
computed from Oa; (c) quadrupole state Uo computed from 
^0 = -Po exp(— r^/2)/27r; (d) quadrupole state U4 after four 
iterations computed from 64,. The insets show the maximum 
intensity Imax = raa,Xx,y \4>\'^ upon propagation, employing 
the approximate solution as an initial condition. 




FIG. 9: (a) 60 — Po exp{—r^ / 2)/ 2tv used to compute the ap- 
proximate dipole- and quadrupole-states in Fig.l^Ja) andlsTc): 
(b) used to compute the approximate dipole-state in Fig. 
ISfb); (c) 04 used to compute the approximate quadrupole- 
state in Fig.|HId). 



(see 0, 03 for details). Then, using for in- 
stance, 128 X 128 mesh it takes only a few seconds 
to obtain the desired solutions. Even if the condition 
6 = Po exp(— r^/2)/27r is not exactly fulfilled, the ob- 
tained eigenvectors may be useful as approximate solu- 
tions. Moreover, it turns out that the following iterative 
process can be used to significantly improve the solu- 



tion. First, from the highly nonlocal approximation of 
the desired solution Uq we compute improved form of 
the nonlocal term 



27r 



(Al) 



and use it to get a new eigenvalue problem Aiv = Xv. By 
solving this new problem, we get Ui which, again, can be 
used in the next iteration, and so on. There is no guar- 
antee that Un converges to the exact solution U, and in 
general it does not. However, by computing the residuum 
of Eq. with U — f/„ it is possible to monitor progress 
of the iteration process. Fig. |S1 shows two examples of 
this technique, performed on the dipole- and quadrupole 
state at power Pq — 1000. The actual 9 used are shown in 
Fig. El In both cases, the initial guess Uq [Figs.|S{a) and 
IHIc)] computed from 60 — Pq exp(— r^/2)/27r [see Fig. 
|5fa)] can be improved by 6 and 4 iterations, respectively 
[Fig-EIb) and|SJd)]. The actual test of the approximate 
solutions is provided by employing them as initial data 
and propagate over a certain z distance. The results are 
presented in the respective insets of Fig. |S1 which show 
that soliton solutions are stable (at least over a propaga- 
tion distance of 2; = 10). Of course, since we treat a non- 
integrable system these solitons are expected to possess 
internal modes [i^ . If we use an approximate solution as 
initial data, some of these internal modes become excited 
and cause oscillations of soliton amplitude. The better 
the approximation the weaker the oscillations. This fact 
is clearly visible when comparing the plots in Fig. |Hta) 
and|HIc), andlHIb) and|SId) respectively. 

These two examples show that it is possible to com- 
pute useful approximate solutions using the above itera- 
tion technique. The big advantage is its efficiency since 
the computations involved take only a couple of seconds 
on up-to-date computers. Moreover, it turns out that 
the technique can be used with non-Gaussian response 
functions, like Eqs. Q and ||SJ), as well. In fact, all solu- 
tions presented in this paper were computed as described 
above. 
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